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^^ ■ Abstract. This paper provides the first provable 0{N log N) algorithms for the linear system aris- 

ing from the direct finite element discretization of the fourth-order equation with different boundary 

^ C^ ' conditions on unstructured grids of size A' on an arbitrary polygoanl domain. Several precondi- 

tioners are presented, and the conjugate gradient methods applied with these preconditioners are 

CO ■ proven to converge uniformly with respect to the size of the preconditioned linear system. One 

main ingredient of the optimal preconditioners is a mixed-form discretization of the fourth-order 
problem. Such a mixed-form discretization leads to a non-desirable — either non-optimal or non- 
convergent — approximation of the original solution, but it provides optimal preconditioners for the 
direct finite element problem. It is further shown that the implementation of the preconditioners can 
be reduced to the solution of several discrete Poisson equations. Therefore, any existing optimal 
or nearly optimal solver, such as geometric or algebraic multigrid methods, for Poisson equations 
would lead to a nearly optimal solver for the discrete fourth-order system. A number of nonstandard 
K^ ■ Sobolev spaces and their discretizations defined on the boundary of polygonal domains are carefully 

^■f-N ■ Studied and used for the analysis of those preconditioners. 
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1. Introduction 

In this paper, we study numerical methods for solving finite element systems for boundary 
value problems for fourth-order elliptic partial differential equations on unstructured grids of an 
arbitrary polygonal domain. The fourth-order partial differential equation has applications in, for 
example, solid and fluid mechanics and material sciences. Many different finite element methods, 
including conforming and nonconforming, have been developed in the literature for discretizing 
the boundary value problems. All these discretizations lead to very ill-conditioned linear systems 
with condition numbers of order 0(h~'^) that are difficult to solve. 

Multigrid methods are among the most efficient techniques for solving these systems. In partic- 
ular, geometric multigrid methods based on a nested sequence of multilevel geometric grids have 
been extensively studied in the literature for fourth-order problems, c.f. [l4ll7ll8l lT8ll28ll331l451l46ll48l 
and references therein. The efficiency of these methods, however, depends crucially on appropriate 
underlying multilevel structures. Because such multilevel structures are not naturally available in 
most unstructured grids in practice, multigrid methods of this type are generally quite difficult to 
use. More user-friendly methods such as algebraic multigrid methods (that can be applied to un- 
structured grids) have also been studied in the literature, see Il5ll6ll23ll30ll34l. But the efficiency of 
these methods applied on the fourth-order finite element problem is limited. Further, there is still 
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no theory to support the methods of these types. In any event, to the authors' knowledge, no math- 
ematically optimal solvers for fourth-order finite element problems discretized on unstructured 
grids are presented in the literature. 

In this paper, we develop a class of methods that fall in between geometric and algebraic multi- 
grid methods. We use the final geometric grid only which has no hierarchy and can be entirely 
unstructured. We circumvent the obstacles inhering in a lack of natural hierachical structure by 
exploring the deep virtue of multilevel methodology in the framework of the Fast Auxiliary Space 
Preconditiong (FASP) method that exists in the literature [|40ll4T| . These methods combine the 
practical advantage of algebraic multigrid methods in that they are easy to use with the solid the- 
oretical foundation associated with geometric multigrid methods. 

To cope with the conforming and nonconforming finite element discretizations for the fourth- 
order problem with different boundary conditions, we present several preconditioners and then use 
the conjugate gradient method on the preconditioned linear systems. One main ingredient in the 
construction of these preconditioners is an auxiliary discretization of a mixed form of the fourth- 
order problem whereby piecewise linear finite elements are used for both the original variable and 
the auxiliary variable. Such simple mixed finite element discretizations often lead to non-desirable 
(either non-optimal or non-convergent) approximations of the original solution; however, a proper 
combination of the solution and some elementary point relaxation iterative methods for the original 
system, such as the Jacobi and symmetric Gauss-Seidel methods, produces a preconditioner that 
can capture the spectrum of the original linear systems well. And, as a result, the conjugate 
gradient method applied to the preconditoned system converges uniformly with respect to the size 
of the systems. The solution of the linear mixed system will be reduced to the solution of several 
discrete Poisson equations. Thus, any existing optimal or nearly optimal Poisson solver, such as 
geometric or algebraic multigrid methods, would lead to a nearly optimal solver for the discrete 
fourth-order system. These solvers work for a relatively large class of finite element systems. 

In order to analyze the optimality and complexity of the preconditioners, a number of technical 
results associated with the Sobolev space H^ n H^^iO.) on the polygonal domain are developed 
in the paper. That is, the finite element discretizations and the continuous and discrete traces on 
the boundary are studied, as the generalization of some existing results for the smooth domain. 
Based on the properties of the newly developed trace spaces and the trace and extension operators 
defined on the trace spaces, the optimality of each solver is shown rigorously on both convex and 
nonconvex domains. 

With various numerical examples also reported to support the theoretical results, this paper 

appears to be the first to present provable 0{N log N) algorithms for the fourth-order equation 
discretized on unstructured grids. 
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The rest of the paper is organized as follows. In the next section, we introduce the model 
problem and basic definitions and notations of finite element methods. In Section [31 we intro- 
duce some preliminaries of the paper, namely a concise introduction to preconditioned conjugate 
gradient (PCG) method and FASP. In Section |4l we discuss several discrete second-order opera- 
tors defined on finite element spaces, and in Section |5l we develop some trace spaces of specific 
Sobolev spaces and discuss their discretization. Sections [6] and |7] are devoted to the optimal solver 
for the finite element problem of the fourth-order problems. Both theoretical results and numerical 
examples are provided to verify the optimality of the solvers. Finally in Section[8l conclusions and 
remarks are given. 

2. Model problems and finite element discretizations 

2.1. Model problems. Let Q c R^ be a bounded polygon, with F = dQ. its boundary, and F, and 
/ = 1 , . . . , ^ as the edges of D, such that F = IJ^i T,. We consider the model biharmonic equation 

(1) A^u = f 

on Q equipped with boundary value conditions of the first and second kind, and its variational 
formulation. In doing so, we find u e M^, such that 

(2) {V^u:V\) = (f,v) VveM,, 

where M^ is a Hilbert space with a certain boundary value condition; namely. Mi = H^{^) when 
the Dirichlet boundary condition of the first kind is considered, and M2 = H^{Q.) fl H^{Q.), when 
the Dirichlet boundary condition of the second kind is considered. In terms of elasticity, the 
first is in accordance with the clamped plate and the second is in accordance with the simply 
supported plate, and they will be referred to thereafter as the first and second biharmonic problem, 

respectively. It is known for w 6 H^(Q.) n //gCQ) that ||Aw||o,n = l|V^w||o,n = \w\2,n = l|w||2,n( GH). 
The well-posedness of the boundary value problem is obvious by the Lax-Milgram Lemma. An 
equivalent weak form of the boundary value problem is to seek u 6 M^, such that 

(3) (Am,Av) = (/,v) VveM,. 

2.2. Finite element problems for model problems. 

2.2.1. Triangulation and finite element. Let Th be a quasi-uniform triangular triangulation of do- 
main Q, Q = UTeTi.T. We denote hr as the mesh size of T, and h as the meshsize. 

For a given triangulation, the finite element space is defined by the elementwise shape function 
space and the continuity of the nodal parameters. In this paper, we study the finite elements whose 

nodal parameters are of the type Na((p) = ■£, (V''"^)(ti, . . . , tkj, where D^ is the integral domain 
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with respect to a, and ti, . . . , t,t^ are k^ identical or different unit vectors to denote the direction of 

the derivative, ka=0, 1,2, In general. Da is a subsimplex of the triangulation; that is, it can be 

a point, an edge, or a triangle. When D^ is a point, the average of integration is reduced to the 
evaluation on the point. Define the degree of the nodal parameter by deg(Na) = deg(a) := k^. 
When deg{a) = 1, there is only one direction involved for the derivative, and the direction is 
denoted by Iq,. 

Let ^ be a subsimplex of the triangulation. Define the neighboring patch of K as cok = UI^"' ^ 

Th : T' Ci K i^ Q}. Given T e Tu, let Pk{T) denote the polynomials on T with a degree not higher 
than k, and Pk{Th) the set of piecewise polynomials that belong to PuiT) at each element T. 

For a given finite element space, the nodal basis is defined to be the set of all the dual basis 
functions with respect to the nodal parameters. Let ipa be the nodal basis function with respect 

to Na- Define coa := [MT : T f] supp(tpa) ^ 0}, #cOa '■= #{T : T f] supp((pa) 4^ 0}, and h^ = 
maxrcw^ hj. Any <x>a can be written as oJq, = Ti |J • • • IJ T^^^^, such that T ^ and T^+i share a common 
edge. 

We make use of the following notations. Let Nh denote the set of all the vertices. A//, = N\ U 
A^^, with N]^ and N^^ consisting of the interior vertices and the boundary vertices, respectively. 
Similarly, let fi/, = &\ IJ &\ denote the set of all the edges, with &\ and &\ consisting of the interior 
edges and boundary edges, respectively. For an edge e, n^ is a unit vector normal to e. For e e &\, 
Tj; and Tf are the two adjacent elements that share e as the common edge, and n^ and nf denote 
the unit outer normal vectors of T^ and Tf on e. Define on the edge e 



dif/ 
5ii„ 



dil/\jL dif/\jR diJ/\jL di//\jR 

r- + ^ = -O^e ■ nt:) + -{n^ ■ n!;), 

dn^ dn^ dn^ dUg 

when e e &[ and e = T^r^ Tf; 

, when e 6 £?. 

I dUe 



2.2.2. Finite element problems. We assume that the finite element space under consideration sat- 
isfies the following conditions: 

Al: the quadratic polynomials are contained in the elementwise shape function space; 

A2: each w/, e M/, is weakly continuous and Vw/j is weakly continuous on each e 6 fi^ in the 
following sense, where n is the normal vector of e and r is the tangential normal vector of 
e: w/, is continuous on at least one point on e, 5„w/, is continuous on at least one point on 
e, and djWh is continuous on at least one point on e; 
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A3: for w/, 6 Mf, and T eTh, 

(4) ||w;,||^^, = Y, hl^''''"'^Kiw,)\ 

Here and after, we make use of <, > and = to denote <, > and = up to a constant. The hidden 
constants depend on the domain. And, when the triangulation is involved, they also depend on the 
shape-regularity of the triangulation, but they do not depend on h or any other mesh parameter. 

Remark 1. The Assumptions Al -A3 are mild, and they hold for most finite elements for fourth- 
order problems, including the Morley element, the modified-Zienkiewicz element /|56|/. the Nilssen- 
Tai-Winther element l^6\l. the Morley-Zienkiewicz. element l^JTil. the Bell element, and the Argyris 
element. 

We define a piecewise H^ function space as 

H\Th) := {w 6 lH^) : w\t e H\T), V T e Tu). 

The operator V^ is defined for w e H^{Th) by (V^w)!^ = V^{w\t). By the assumptions Al and A2, 
the functional | • b,/, defined as |w/,|2,a = l|V/^w/,||o,n is a nontrivial semi-norm on M/,. 
In practice, the finite element space M/, i that is associated with H^(Q.) is defined by 

Mh,i = Iwh 6 Mh : Naiwh) = 0, when deg(a) < 1 and Da c dOl. 

And the finite element space M/,_2 that is associated with H^{Q.) n H^iO.) is defined by 

Mh,2 := Iwa 6 M/, : -f {drf'^^^^Wh = when D^ c dQ, and deg(a) = 0, 1 

By definition, M/j i c M/,_2. The semi-norm | ■ I2,/, is a norm on M/, 2- 

The finite element discretization of © is to find ui, 6 Mh^t {k = 1,2) such that 



I fVh, 

Ja 



(5) ah(uh,Vh) := (V;,m/, : V,;vft) = | fvh, Vv/, e Mia- 

We introduce the operator A/, ^t : M/,_;t ^ M/, ^t, 

(Ah^kVh, Wh) = ah{vh, Wh), V Vh, Wh e M/,,^. 
Then the linear system to solve is Ai^j^ui, = fh, with //, the L^ projection of / on Mhi,. 
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3. Preliminaries: PCG and FAS? 

3.1. Preconditioned conjugate gradient (PCG) method. The PCG method is the basis of all 
the preconditioning techniques studied in this paper. For a linear system Au = f in which A is a 
symmetric positive definite (SPD) operator, the PCG method can be viewed as a conjugate gradient 
method applied to the preconditioned system BAu = Bf. Here, B is an SPD operator, and BA is 
symmetric with respect to the inner product (•, Ob' := (B^^-, •)• 

Let Uk,k = 0,1,2,- ■ ■ be the solution sequence of the PCG algorithm. It is well known that 

r^^ II II ^./ V^<^-1 \'„ „ 

(6) \\u-Uk\\A <2 — - M|m-mo||a, 

\^/KiBA) + 1/ 

which implies that generally the smaller the condition number k(BA), the faster the PCG method 
will converge. However, though it is sufficient for many applications, the estimate given in Q 
is not sharp. An improved estimate can be obtained in terms of the eigenvalue distribution of 
BA lfTl fT7ll38l . More specifically, we focus on the case in which we can divide cr(BA), the spectrum 
of BA, into two subsets, cro{BA) and (Ti(BA), where the eigenvalues in cri are bounded from above 
and from below and ctq consists of all the "bad" eigenvalues remaining. The lemma below can be 
found in many references, e.g. tlj. 

Lemma 2. Suppose that (t(BA) = crQ(BA) U (Ti(BA) such that there are m elements in cro(5A) and 
A 6 [a, b]for each A e cri{BA). Then 



(7) ||„_„.„,^24^^) 

\ Sbja + 1 / 



k-m 

\\U - MolU 



where 

K= max n I1--I- 

ieo-i(BA) 11 Li 

liea-Q(BA) ^ 

If the eigenvalues of BA are ordered as < Ai(BA) < • • • < An-mi,(BA) «; /l„_„,^+i(5A) < • • • < 
An(BA), then the convergence rate estimate dV]) becomes 



(8) J. ^ < 2 

\\U - UqWa 



( JA.,^„JBA)/AABA) - 1 f """ 

, k > mh. 



Uk\\A ^^ yJAn-m,{BA)IAy{BA) - 1 



\^lA:;ZjBA)iM{BA) + \, 



In this case, /i„-„,,/5A)//li(5A) plays a dominant role in determining the convergence rate, and it is 
sometimes called an eff"ective condition number of BA. In general, let A be an SPD operator, and 
the eigenvalues of A be ordered by < A\{A) < A2{A) < • • • < A„{A). Define the m-th eff"ective 
condition number of A by /cjjf (^) •= ^n-m{A?)l A\{A). A similar consideration can be found in [|43l . 
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Obviously, K^f(A) = k{A) and Kf^(A) decreases as m grows. Then based on ([8]), given a tolerance 

e < 1, the number of iterations of the PCG algorithm needed for ||"~"*H^ < e is 

° ll«-«olU 

(9) k>m„ + log s(K^^(A)y'\ 

To estimate the effective condition numbers, a basic tool is the Courant minimax principle (see 
e-g-5 [Mil), which leads to the following lemma: 

Lemma 3. Let V be a Hilbert space and A : V ^ V be an SPD operator on V. Suppose the 
eigenvalues of A are ordered a* /li < /I2 < ■ • • < /t,„ then for any subspace Vq <zV with dim(VQ) = 
m, the following estimates hold: 

(10) /i„, < max and /l„+i_„,(^) ^ mm . 

veV'o (V, V) veVo (v, v) 

3.2. Theory of the Fast AuxiUary Space Preconditioning (FASP) method. In this section, we 
give a summary of preconditioning techniques based on auxiliary spaces as developed in [|20[|39l - 

m. 

Let V stand for a real Hilbert space with an inner product a(-, •) and a (energy) norm || ■ |U. Write 
A : V ^> V for the isomorphisms associated with a(-, ■)• Let s(-, •) be another inner product on V, 
and 5 : y 1-^ V the isomorphism associated with s(-, •)• Let Wi,. .., Wj, 7 6 N be Hilbert spaces 
endowed with inner products dj(-, ■), j = 1, . . . , 7. Write Aj : Wj 1-^ W- for the isomorphisms 

associated with a/-, ■), j = 1, . . . , 7. Furthermore, for each Wj, we need a linear transfer operator 
Uj : Wj H^ V. We tag the adjoint operators by "*". 

Then the fast auxiliary space preconditioner (see [@0]5D) is defined by 



(11) B = R + Y^njoA]^ oU*j, withR = S 



Theorem 4. l \2(J\\40\l Assume that 

(1) There are constants Cj > 0, such that 

WUjWjWa < CjajiWj,Wjf\ -iWj G Wj-, 

(2) There is a constant c^ > 0, such that 

\\v\\a<cXv,v)"\ VveV; 

(3) There is a constant cq > 0, such that for any v £ V, there are vq g V and Wj £ Wj such that 

J 

V = vo + ^ UjWj, 
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and 

J 
5(vo,Vo) + Y^aj{Wj,Wj) < CqIIvII^. 

7=1 

Then 

k(BA) < cl[cl +cl + --- + cj). 

In applications, V and all Wj are usually finite element spaces with bases consisting of locally 
supported functions. Plugging basis functions into the bilinear forms would lead to the algebraic 
representation of the preconditioner. 

4. Discrete Second-order operators 

The Laplacian operator and the Hessian operator are second- order operators of fundamental 
importance. In this section, we study their discretization on finite element spaces. 

For a linear space M and its subspace M* c M, denote codim{M^„ M) as the codimension of M* 
in M. By definition, we can obtain the following lemma directly. 

Lemma 5. Assume a linear space M and a linear operator L on M. For any subspace M^ of M, 
codim(L(Mt), L(M)) < codim(M^,M); ifker(L) c M*, thencodim(L{M^),L(M)) = codim{M^,M). 

Assume 'R{^) = A{H^(Cl) n H^iQ)) := {Aw: w e H\£l) n //^(Q)}. Following, e.g., U6J, we 

have codim{K{A), L^(Q.)) = mq, where Mq is the number of reentrant comers on F throughout this 
paper, and by a "reentrant comer" we refer to the comer whose interior angle is bigger than n. 
When the domain Q. is convex, /ng = 0. 

4.1. Discrete Laplacians on linear finite element spaces for H^(Q.). Given a triangulation 7~/„ 
denote by V/, c H^{Q.) the continuous linear element space on 7"/,, and V/,o = VhC\Hl^(Q.). Moreover, 
define Vf,b as the complementary subspace of V/jq in Vf,, namely, V/,/, c Vf, and Vf, = V/jo © Vhb, and 
Bf, as the trace space of V/, on F. Define Sf, as consisting of functions that vanish on the comers of 
F, and S^ the complementary of S/, in S/,. 

Definition 6. Define the discrete Laplacian operators of first and second kind on V/,o as 
(12) - Ah,k ■ Vho -^ Via such that {-Ah,kWh, v/,) = (Vw/,, Vv/,), V v/, 6 Vh.k, "^Wh^ Vm; 

where Vh,\ = V/,, and Vh^2 = Vho- 

Define Vho := [ph 6 Vm : AnaPh 6 7?(A)), then A^.z^o = 'R(A) n V^. As ker{Ah,2) c Vho, 
codim{Vko, V/,o) = codimiAh^Vho, Ah^Vho) = codimCRiA) n Vho, Vho) < mo. 
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Each of the two Laplacian operators has a second-order equivalent description in terms of nor- 
mal derivative jumps on element edges. Let £/,j = fi/,, and &h;i = &\. 

Lemma 7. The Laplacian operators can be described by 

(13) IIA,aP/Jlo,n ^ Yj ^^"' /l^^f ~ W^KkPhWln^ 

which holds for all pu e V/io and k = 1,2, with this exception: when k = 2, the right inequality 
above only holds for ph 6 Vm- 

Proof Let ph 6 V/,o and k = 1,2, then 



e^£>h,k 



Namely, ||A,,,p,||g„ < ^ h~' J \-^\ . 

Now we turn to the right inequality of (IT3T ). We first consider the case k = 2. Let p e H^{0.) be 
the unique solution of 

-Ap = -Ah,2Ph, infl, 
/> = 0, on 5^. 

Then, when Ah^2Ph e "^(A), it holds that \ph - p|i,n < hWph.n < h\\AhaPh\\o,ii- Thus, 

z *.- filer =?*--Ii^r 

^14^ ^E J](hf\p-p,.\lT + \p-Pn\lT) 

ee&,,2 dTDF 

< \p\la < l|A,,2;'Allo,n- 
The right inequality of (fT3]) is proved for k = 2, and we turn to the case ^ = 1 . As 

l|A/,,ij!?/,||o,n = sup — - — = sup 



qi.eViMO] \\lh\\o,n qhSViMO] WQhWo, 



a 
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we will find some r/, 6 V/, \ {0}, such that 



^"'IM) 



2 ^ (V p/„Vr/,)' 



< 



If Q is convex, then by the same argument as for fc = 2 with the homogeneous Dirichlet problem 
replaced by the homogeneous Neumann problem, we obtain that \ h~^ \\~n~\ ~ il'^/iJ/'/illon- 

Then (fTSi) holds for r/, := -A/, ip/,. 

If ^ is not convex, let Q. be the convex hull of all the comer points of Q, and let 7/, be a quasi- 
uniform triangulation on Cl such that 7"/, is a subtriangulation of T/,. Similarly, the set of edges 
£ft is defined on 7"/,, and the discrete Laplace operator A;, j is defined on V/,, the continuous linear 
element space defined on T/,. Denote by ph the extension of ph to V/„ such that Ph\a,\a. = 0- Then 



i^- im -%'- [i^s 



esSh.i eeSi, 



the shape regularity of 7^, only. Define r^ := r/,|n so that 



Denote r/, = -Af^yp^, then / ^^.M "^^ ^ CiHr/JP - for a constant Ci depending on Cl and 

ee£i, 



eeSh 






Therefore, (fTST i holds for this r/,. The right inequality of (IT3]) is proved for fc = 1, and this finishes 
the proof of the lemma. n 
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4.2. Discrete Hessions on finite element spaces for H^(Q.). Let M/, be a finite element space for 
H^{Q) as described in Section 12.2. 2[ Define interpolation operators on finite element spaces as 

h : Mh2 -^ Vho, (IhWh)(a) = -^( V wMa)), Va e Nj^, and Vw/, e Mhx, 

nft,i : Vho -^ Mhj, Na(U,,^\Ph) = 0, when D„ c dQ; 

Na{^h,\Ph) = - — y NaiPhW), else; 



n/,,2 : V/,0 ^ Mft^2, Na{^h,2Ph) = T — ^ T 5n(;'ftlr)(n ■ t„), when D„ c (9Q, and c/eg(Q;) = 1, 
NaiUhaPh) = - — ^ Naiphk), else. 



#^- re. 



Define M/,,2 := {w/, 6 M/,,2 : Ahahwh 6 •R(A)). Then, wi, 6 M/,,2 if and only if //,w/, G Vho- By 
Lemma[5l codim(Mh2,Mh2) = codim(IhMh2, hiMij2) < co(i/m(A;'^2^/jo, V/io) < wq. 
Theorem 8. We have a stable decomposition ofMi,/,: 

(16) WlwhWln < Yj ^T^W^h - Tlh,khwh\\lj + i^hJirWH, Ah,khwh) < WlwhWl^, 

Ten 

which holds for all W;, £ Mk,k cind k = 1,2, with this exception: when k = 2, the left inequality 

above only holds for Wh e M/,_2. 

Proof The theorem follows from Lemma 171 and Lemma [T2] below. n 

Lemma 9. Let m ^ I be an integer The following equivalence for y,ySi, . . . ,/3„, 6 R depends on 
m only: 

.. m .. m 

- Yiy-Pi? = (r -fif + - Yifii - A+i)', 

m ^—^ m ^—^ 

(=1 (=1 

where fi = ^{JSy + ■■■ +J3,„), A„+i = ySj. 

The proof of Lemma |9] is straightforward; therefore it is omitted. 
Lemma 10. / I?7l/ It holds for all Wh 6 M/,,2 that 

(17) IK - hWhWlj ^ Yj ^T'Mlr, Vr e T,,. 

T'dWT 

Lemma 11. For all ph e V/,o, it holds fork =1,2 that 

(18) Yu'^r'm^'^PH-p,.\\lr<Y,K' fll^f- 

le/i, eeiju.k 
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Proof. By Assumptions Al and A3, we have 
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(19) 



TeTh deg(a)=l Tccj„ 



Let Na be a nodal parameter with deg{a) = 1, and we rewrite the element patch oJq, as cOa = 
{^n!=i"' such that T" and T" share a common edge, j = 1, . . . , #a)a - 1- Denote r#^„+i = Ti, then 

(20) Yj \^-iPh\T) - Na{UH,kPk)f 

1 |2 *'^" 2 

K(nh,kPh) - ^ — ^ K,iPh\T)\ + ^ (A/'„(pftlrj) - J^a(Ph\rjj) , k=l,2. 



TCOJa 



= #COr 



#a)a 



TCdla 



7=1 



If Dq, <t dCl, then by the definitions of Tih,k, 



#(i)a 



(21) 



2] \!^a{Ph\T) - Na{Tih,kPh) = ^ [NaiPhlrj) - NaiPhlrjj) 

;=i 



TCdJa 



and by the continuity of the piecewise linear function ph, we obtain further that 
(22) J] |A^.(PAlr) - A^a(n,,,;?,)|' < _^ ^^.^' fl^f ' k=\,2. 



Tca)„ 



enDai^Q 



If Z)q, c dCl, then for 11/, i we obtain that 



(23) 



J]\KriPl,\T)-K(n,,,p,f = J]\N,ip,\Tf < Yj K' {{^t- 



e-drie 



dPhi^ 

\J^aKPh\T)\ S: l_^ '^e \ I 

For H/, 2, there are two cases. If deg(a) = 1, then noting that p/,|r = 0, we obtain 

^ \Na(Ph\T) -Na{TVh,2Ph) 

= Yj\\ ^n(P/'lr)(n • tJ + r dr{Ph\T){T-ia)-Na{Ilh,2Ph) 



TEOJry 



TCOJa 



Da 



(24) 



< 



Elf «" 



(Pft|T)(n • t„) - NaiUhaPh) 



= £( f d„(pH\T,) - f d„ip„\T^jf(n ■ tS- + 






(P/ilr)(T • t„) 



(;?/ilr)(T • t„) 



enAr*0.e(tf5n 
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Otherwise, by definition, we have 

(25) Y, \^»(Ph\T) - NaiUh,2Ph) = Y, (NaiPklrj) - NaiPkkjj) . 

Tca/„ 7=1 

And, by the continuity of ph again, we obtain that 

(26) Yj \^-(P'r\T) - K.(U„2Pk)\' < Y ^'^' jl^f- 

Summing all the inequalities above leads to (fTS]) fork= 1,2. 
Lemma 12. It holds for k = 1,2 that 

(27) IIV^w.ll^^ = Y ^T^Wwh - ^KkhWhfoj + Y ^'^"' fl^f ' ^^^^'' ^ ^'''^- 

Ten ee£,,, ^' ^^^ 

Proof. Denote w]^ := hwh for any Wh 6 Mh,k, k = 1, 2. By inverse inequality, we have that 

Ten Ten Ten 

Thus, by Lemma [TOl and Assumption A3, 

Ten 



(28) ^ E ^^^' Z /^?'''^'^'^tA^«(>^/>) - NM\t)) 

Ten (^a^T 

= ^ /i2^^^(-)-2 ^ (yv„(w,) - yv„(w;jr))'. 



TCOJrr 



It is straightforward to obtain that, for any a, any wi, e M^^t and any p^ 6 V^o, ^ = L 2, 

(29) ^ (W,(w,) - K(ph\T)f = #aj,(K(Wh) - NaiUH,tPh)f + Y (^«(n/ai?ft) - KiPh)f . 

Then by Assumption A3 again, given any Wh £ M/^,^, 



Tc(j)„ 



(30) +(A^.(n,X)-^«(<lT))] 



Ten Ten 



\wh\lh 
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therefore, by (fTST ). 

'■." 1 1 

On the other hand, as 
we obtain 



E ^= I[0 * E ^^- 1[^4^]^ ? ^- Ii^r ^ «i.. 



2 ^r'llvv. - n„kwi\\lr + J] K' I [^] , fc = 1, 2. 



Combining these two points, we obtain (|27l) . n 

5. Trace spaces on the boundary and their discretizations 

5.1 . Trace spaces related to H\Q). As the trace space of H\Q.), H^'^{Y) is a Hilbert space with 
respect to the norm ||/l||i/7r = inf llwllin. Denote H'^^^{Y) as the dual space of H^^^(Y). 

weHH^),A=w\r 

Given A e H^^^(r), there exists a unique ua 6 H^(Q.) satisfying 

-Auji = in Q 
ua = A on dQ. 

This defines a harmonic extension operator by EA = u^. Then, H^/llli^n = iUlli/2,r for A e H^^^iY). 

For any edge F, of Q, define the space HH^iYi) := {/l 6 L^(F,) : 1 e //^^^(F)} where A is the zero 

extension of/l into F\F,-, and the space is a Hilbert space with a norm given by ||/l||^i/2,px := ||l||i/2,r- 

Moreover, H^^iYi) is the interpolated space halfway between the H^iYi) and L^(F;) spaces. We 



refer to [1221 and references therein for more details. 

5.2. Trace spaces related to H^{Q) n Hq(Q.). We consider the following subspace of H^'^(Y): 

(31) Hl'\Y) := {A 6 H"HY) : A^i e //^^'(F), / = 1,2, . . .,K}, 

where xt is the characteristic function on F,, with the norm 



cX-=(J]UXi\\l2,ry''- 



11/2. 

!=1 
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Proposition 13. For any A e ^^(F), the sufficient and necessary condition of A e H^ (F) is 

A\r, 6 H]l^{Yi)for i = 1,...,K. Moreover, ||i||i/2,c,r = ZL UlrMH'i^r.y ^^^ Hl'^iY) is a Hilbert 

space with respect to \\ ■ ||i/2,c,r- 

Theorem 14. The following identities hold both algebraically and topologically: 

(32) Hl'^iY) =\a£ L^CF) : An £ {H^I^{Y)f I = U 6 L^{Y) : At e {H^'^{Y)f \. 



More specifically, IUIIi/2,c,r = IUn||i/2,r = \\AT\\ii2xfor A e Hl'^iT). 

Proof. To begin, we prove the first identity in (|32] |. Let A e H^^iF). We have An = 2^i ^^Xi^Tr 
By definition, Axt 6 H^'^iY); hence, AxiUr, e (H^'^(Y)f; thus, ^n 6 (H^'^(Y)f. 

Now let A G L^(F) be such that An e (//'''^(F))^. It can be easily verified for j = i - 1 and / + 1 
that Axi = An ■ ^ ^ on F,- U Yj and Axi = An ■ ^ ^ = on Fy. Namely, Axi can be expressed as 

' j i J 

a linear combination of function and functions An ■ _\ for j = z + 1, which belong to H^^^(Y), 

I,- Tlj 

with the coefficients in C°°(F). Therefore, Axi itself belongs to H^^^(Y). 

It is easy to verify that ||/ln||i/2,r is a norm on /I e //^ (F). Moreover, ||/ln||i/2,r ^ IUIIi/2,c,r- It is 

1 /7 

straightforward to verify that H^. (F) is a Hilbert space with respect to both the norms || • n||i/2,r 

and II • ||i/2,c,r' therefore, by open mapping theorem, ||/l||i/2,c,r = IUn||i/2,r for A £ hI^^{Y). This 
finishes the proof of the first identity. 

The second identity in (|32)) is obvious given that t = n^. And, the theorem is thus proved. n 

Denote H'^^^iX) as the dual space of //^^(F), with the norm ||;tll-i/2cr = sup — — ^ , 

where <•, ■) is the duaUty between H~^'^{Y) and hI''^{Y). 

Define a space of biharmonic functions H^(Q.) := {u £ H^(Q.) n H^(Q.) : (Au,Av) = 0, Vv G 
ifo(Q)}. It can be verified that V(H^(Q.)) is a Hilbert subspace of H(div;Q.) and A(H^iQ.)) is a 
Hilbert subspace of L^(Q). In the remainder of this section, we show that hI'^(Y) and //^ (F) are 
isomorphic trace spaces of V(//^(Q)) and A{H^{£l)), respectively. 

Theorem 15. H^iY) is the normal derivative trace space ofH^{Q.) n H^^iO.) in the sense that: 

(1) Ifu £ H\a) n //i(Q), then ^^ £ H'J\Y) and |||||i/2,c,r < ll^l^.n. 

(2) Given any A £ hJ (F), there exists a unique u £ H^(Q.), such that |^|r = A and ||M||2,n ^ 

i/2,c,r- 



OPTIMAL SOLVERS FOR FOURTH-ORDER PDES 17 

Proof. Given u 6 H^'D.) n //^(Q), then Vm e (//H^)) and V^lr e (//^/^(r)) . As f^ = along Y, 

it holds that |n = Vu and | e Z/^'Cr). Moreover, |||i||j/2,,,r = ll^"lli/2.r ^ Il"ll2,n- This proves 
the first part of the theorem. 

To prove the second part, we consider an auxiliary Stokes problem for any A e H^^ir): 

' Ai/f + Vp = in Q, 

V • i/r = in Q, 

i// = At on r. 

This problem, thanks to the trivial fact of r • n = on F, obviously admits a unique solution, 
such that (ij/,p) 6 (h\Q.)) x LI(Q) satisfies Hi/'lli.n ^ IUT||i/2,r- (^^^ Theorem 5.1 and Remark 
5.3 on pp. 80-83 in [fT2l|'). As V • i/r = 0, there exists awe H^(Q.) such that i/r = cmlu (see 
Theorem 3.1 on p. 37 in [HIl). And, |^ = curk ■ t = iff ■ t = A on Y. Given v G HgCfl), 
(Am, Av) = (V X curlw, V x curlv) = (V x i/r, V x curlv) = (V x i/r, V x curlv) + (p, V • curlv) = 0. 
Further, |^ = VM-T = i/f-n = along F; therefore, m is a constant along F, and thus we may choose 

u £ i/o(Q). Then, u 6 H^(Q.) and ||M||2,n < Iwh.n + \u\2,d. = l|VM||i,n = l|iAlli,n ^ IUT||i/2,r = l|/l|li/2,c,r- 
The uniqueness of such a m is straightforward. This finishes the proof of the theorem. n 

By Theorem [HI given A 6 //^^(F), there exists a unique ua 6 H^{Q) such that ^ = A. This 
defines an extension operator E^ : H, (Y) -^ V(H-(Q.)) by E^/l := Vm,i. Moreover, we can define 
a trace operator T^ : ViH^iQ)) -^ hI'^(Y) by Tf{Vu) = Vw ■ n for w e fl^2(^) jhe following 
lemma follows from Theorem [T?] directly. 

Lemma 16. T,"' /i' an isomorphism from (V(H^(Q.)), \\ ■ ||,/,v) onto hJ (F), E^ is an isomorphism 
from hI'^{Y) onto (V//^(Q), || • \\dh), and Tf o Ed = Id^^m^^^ and E^ o Tf. = Id^^i^^))- 

Remark 17. It is easy to verify that for u 6 H^{Q.) n H^{Q), \Vu\div = |VM|i,n, therefore, || • ||i,n and 
II • \U. are equivalent on V{H\Q.) n H^iO)). Let ff : ViH^{Q.) n //^(Q)) ^ H^c'^(Y) be defined by 
Tf(Vu) = Vu ■ nfor u 6 H^(Q.) n H^iO). Then the range ofTf is hJ (F), and Tf is a continuous 
extension of Tf. 

Theorem 18. The spaces H~ (F) and A(H^{Q.)) are isomorphic in the sense that 
(1) For any w G H^(Q.), there exists a unique ^„, e He (F), such that 

(33) {U ^) = (Aw, Av), V V 6 H\a) n //^(Q), 

an 
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and ||^wll-i/2,c,r ^ l|Aw||o,n. 
(2) For any ^ e H^ (0- there exists a unique Wi; 6 H^(Q.), such that 

(34) (Aw^, Av) = <^, ^), V V e H\a) n hI{^), 

anJ||Awf||o,n < ||^||-i/2,c,r- 

Proof. Let us first prove the first part of the theorem. Given A e He (F), then by Theorem 
there exists a v,; G H^{Q), such that A = ^Ir- We define a Unear functional ^„ on H^. (Y) by 
{^,,,A) = (Aw, Av,i). For any fi e hI'^{Y), denote v^ e //^(Q) n HliO.) such that 5^|r = // and 
l|V;j||2,n < C||/z||i/2,f,r with a generic constant C. Then, 

iirii (^.)") (Aw^,Av^) 

ll^ll-i/2,.,r = sup -— = sup —— < C||Awf||o,n. 

To prove the second part of the theorem, let ^ 6 //^ (F). Then, there is a unique w^ 6 H^(Q) 
such that (Aw^,Av) = <^, |^), Vv e H\Q.) n //^(Q). Furthermore, (Aw^,Awf) = U,^) < 

ll^ll-i/2,f,rll^lli/2,.-,r < ll^ll-i/2,c,rl|Wf||2,n < ll^ll-i/2,f,rllAwf||o,n; therefore, ||Awf||o,n < ll^ll-i/2,c,r- 
The proof is finished. n 

By Theorem [TSl given ^ e H^ (F), there exists a unique w^ e H^(Q) such that (Aw^, Av) = 

<^, ^) for all V 6 H^(Q.) n //^(Q). We define an extension operator E, : H~^^^(Y) -^ A{H^{Q)) by 

E,^ := Awf. Moerover, we define a trace operator T,^ : A(H^-iQ.)) -^ H~^^^{Y) by <r,':(Aw), |^) = 

(Aw, Av) V V e H^iQ) n //^(II) for w e ^^(Q). The following lemma follows from Theorem Ml 
directly. 

Lemma 19. T^' i5 an isomorphism from (A(H^(Q.)), \\ ■ ||o,n) onto H^ (F), E^ ^■s' an isomorphism 



/rom He (F) onto (A(// (Q)), || • ||on), and T'^ o E^ = W^-i/2,pj and E^ o T^ = Id^^^fj: 



(m- 



It is easy to verify that if i 6 H^'^{Y), then {EA, Av) = (E,A, Av), V v e H^Q.) n //^(fl); that 
is, EcA is indeed the L^-projection of EA in 'R(A). In particular, when "^(A) = L^(Q.), which 
holds when Q is convex, E^ and E coincide on H^^^(Y). Therefore, E^ : Hc^'^(Y) -^ !^(A) is a 
generalization of the harmonic extension operator to a larger function space H~ (Y) D H^^^(Y), 
and £■,. is called a generalized harmonic extension operator. 
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5.3. Generalized Poincare-Steklov operator and its inverse. Given A e H],'^{Y), let w,i G 
H^(Q.) be such that Vw,} = E^jA. Define a generalized Poincare-Steklov operator as T^ : hI (T) -^ 

H-"\T) by TA := r,f(Aw,). Then 

(35) {T,A,n) = {VE,iA,VE,^), \/ A,fi e hI'^Y). 

Given ^ e H~^'^{Y), let w^ e H^{Q.) be such that Aw^ = E^^. Define an operator 5, on H;}''^{Y) 
hy Sc/^:= Tf{Vw^). We directly obtain that 

(36) {^,SS = {E,^,ES, 'i^,^eH-'l\Y). 

Lemma 20. The operators Tc and S ^ are both algebraic and topological isomorphisms between 
He {Y)andHj (F), andT^oSc = Id^-nify^and S c°Tc = Id^in,^.. Moreover, {^,ScO ~ ll^ll-i/2cr 
for ^ 6 H;"\Y) and {T,A, A) = UW^^^^^^for A e hI'\Y). 

Proof. Firstly we show that T^. is a bijection. It is straightforward to verify that T^ is injective. Now 
let// 6 H~^'^(Y), then there exists a w^ 6 H^iQ), such that Aw^ = Ecfi. Denote A = TfVw^, then 
by definition yu = T^A. Thus T^ is surjective. Similarly, 5^ is a bijection. By definition and the fact 
that Ec and T^ are inverse to each other and Ed and Tf are inverse to each other, TcoSc = Id^-\i2^^. 
and 5 c oTc = Id^iii,yy 

For any ^ e H~^''^{Y), it can be derived by Theorems [B] and [H] that \\S ,^\\ii2,,j- < ||^IUi/2,c,r- 

Similarly, we can prove ||rf/l||-i/2,f,r ^ imii/2,c,r for A e H] (Y). 

The proof is finally finished by noting (|35T ) together with Lemma [19] and ( |36l ) together with 
Lemma [T6l n 

By means of the inverse generalized Poincare-Steklov operator 5 c, the first biharmonic problem 
can be decomposed to two second biharmonic problems, as in the proposition below. 

Proposition 21. Let u e H^{Q.) solve (Aw, Av) = (/, v) Vv e H^i^i), then u can be obtained by 
seeking {u, ^, u) 6 {H'^{Q.) n //^(H)) x H'.^'^iY) x (//^(Q) n //^(Q)j, such that 

(37) (AS, A;.) = (/, p), \/pe H\Q.) n hI{Q.), 

(38) <^,s,r) = -(As,£,r), Vr e //^^/'(r), 

(39) (Am, Av) = (Aii + £,^, Av), Vv 6 //^(ri) n //^(II). 

Proof Let y e //;'^^(F), then by ^ and dMI), we have (Aw, E^y) = (M + E,^, E,y) = 0. Namely, 
<7, |) = (E,y,Au) = for any y 6 //;^^^(F), thus f^|r = 0, and u 6 HliO). Further, for any 
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V e HliQ.), 



(Am, Av) = {M, Av) + {E,^, Av) = (/, v) + <^, ^) = (/, v). 

an 



This finishes the proof. n 

5.4. Linear element spaces on the boundary. Given a subset F (Z Q., define a restriction opera- 
tor as Ip : Sh ^> Sh by IpV{x) = v{x) if x is a vertex in F, and IpV(x) = if x is a vertex out of F. 
The lemma below is about the stability of I p. 

Lemma 22. /I?4l/ Let F be a vertex on dQ or an edge ofdQ.. then ||/°/l/,||i/2,r ^ (1 + 1 log/2|)||/l/,||i/2,r- 
Lemma 23. It holds for Ah g S/, that 

/r\ II 1 II {^h,fJ-h) 

(1) Uh\\-\!2.T = sup —- ; 

IJi,eSk \\Mh\\l/2,r 

(2) |U„||o,r < h~''^\\Ai,\U/2x; 

(3) p;il|-i/2x</^'^'(l + |log/i|)|l4lkr- 

Proof. Let Qh be the L^ projection to S/„ then ||<3/,/l||i/2x ^ Plli/2,r for A e H^'^iY). Therefore 
ll^ll-i/2,r = sup,,^..(r) Ml < sup^,^..(,) i^ < sup,,.,^, ^^ < |U.|i_i/2,r- This proves the 
first item. The second item follows from the first item and the inverse inequality. 

Now, any ///, e S,, can be decomposed as yU/, = fin + jdl- As IKIIi/2,r > IIAJIi/2,r - IKIIi/2,r 
and ||ju^,||i/2,r < (1 + |log/z|)||yU/,||i/2,r, we obtain that ||y&/,||i/2,r + IK;ili/2,r < (1 + |log/z|)||jU/,||i/2,r. 
Then, as fi,, vanishes at the comer, (4,y&/,)r < Ul\\o,r(h\\fih\\i,.n,pp(Ai)) < /?^''^IU/'llo,rlly"/illi/2,r. and 
(^IfiDr < ll4llo,r(/^IKIIi,r) ^ /^'^'ll4Jlo,rlKlli/2,r, we obtain that 

P^ll-i/2,r = sup -^ < (1 +\logh\) sup -^ „;„ < h"\l + |log/i|)||4l|o,r- 

p/,eS/, ll/^Alll/2,r iii,eSi, ll/^/,lll/2,r + llA'/illi/2,r 

The last item is proved, and the proof of the lemma is finished. n 

Lemma 24. Let F,- be an edge ofdQ. Denote S/,,- := S/Jp, = ^k n //o(r,), and define Qm the I? 
projection to S^- Then, 

(40) lU,- - a.-^,-||o,r, + h"^\\QHiAi\\^^^2^^_^ < h^^UtW^m^^^, V A^ 6 Hl^^Yd. 

Proof. Firstly, it can be proved that 

Ui - QkiAiWox^ + m,nAi\\o,r. < IU,llo,r„ V i,- 6 L'(F,-); 



^^^^ ' Ui - QhMox. + /^lia«vl;||i,r, < h\\Ai\h,r., Vi,- 6 //^(F,). 
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Then, as H^^iYi) is the interpolated space halfway between the Hq(Yj) and L^iXd spaces, the result 
is obtained by interpolation. n 

Lemma 25. Define Qh as the t} projection to S/,. Then 

(42) U - QhAWv + h'l^\\QHAhi2,c,Y < h'l^Uhi2,c,T^ Vie H'J^Y). 

Proof. For any A e H}J^(r), by definition, (QhA)\r = <3ft,(^lr,)- Denote A, = i|r, for / = \,...,K; 

then, |U||i/2,.,r = 2^i Uihl^ir^y lia,^lli/2,cr = Zf=i WQhMH'^n^r.) and \\A - Q,A\\o,t = Zf=i Ui - 
^/«'^(llo,r,- Then summing (ffOt from / = 1 to / = ^ leads to (|^ . This finishes the proof. n 

o o o 

For Ajj e Sh, denote its decomposition by A^ = Aj, + A^, with A^ e Sf, and A^^ e S^, such that 
Ah(x) = if jc is a comer point, and A'f^(x) = if jc is a vertex other than the comers. 

Lemma 26. For any A^ e Si„ \\A,,\\ox < h-^'Y^„\li/2,cT and |l4IUi/2,c,r < h^'^U^lWr- 
Proof. By definition of HJ (F) and inverse inequality, we have 



K 

II 1 ., ,, 

(i\\\l2,T 



2°° ° ° ° \~' ° 

lUftllox = (^ft' '^h)T < m/!ll-l/2,c,rP/illl/2.f,r = P/ill-l/2,f.r ^ \\^hXi\ 

i=\ 
K 

< UhW-m^cxYih-'lY^nXiWox) ^ /^"'^'ll^ll-i/2,c,rlUAllo,r- 
(=1 

Thus, ||A/j||o,r ^ ^^'''^IU/ilUi/2,f,r- To prove the second part, for any ip G //^^(F), denote ^/, = Q^^p 
and(4,(^)r = (i;,,^,)r + (4,^ - ^/,)r < /^'^'ll4llo,rll<^/>lli/2,c.r + /^'^'ll4llo,rlMli/2,.,r. Therefore, 
Il4ll-i/2,c.r = sup^g^i/2(p) ii^jj^ < /z^^^ll^lloj. This finishes the proof. 

Lemma 27. It holds for Ah e S/, that 

\\UUi2,x + h''^\\A,\\o,r = Uh\\-i/2,c,r + h"H\4\\o,r. 
Proof. On one hand, 

,1/2||}.|| „ ^113.11 . ,. ^ I a1/2|| 1.11 L incii . ,_ „ I ;,l/2||ic| 



D 



ll^/,ll-i/2,.,r + h'"U,r\\o.r < P/,II-i/2,.,r + h' ' 'UhWox + Il4ll-i/2,.,r + h^'^UlUr 



<\\M-u2,cx + h''^Ui\\or; 



on the other hand. 



\\h\\-u2,x + h'/^UlWox < IUftlUi/2,.,r + IU;,ll-i/2,c,r + h^'^UlWox 

< Uh\\-U2,cx + h^^UlWox < \\M-i/2,cx + h'^^UnWox- 



The proof is thus finished. 
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Lemma 28. It holds for Ah 6 S^ that 

(43) Uu\Ui2,cx + h'i^Uh\WT ^ Uk\Ui2.T < m\-mx + (.^ + \^ogh\)h'i^Ul\WT. 
Proof. By inverse inequality, /z'''^m/,||o,r ^ lU/JUi/i.r- By definition, 

,... ., ,,, (^,ju)r / (/l/„//)r ^ (4,/^)r i, , i, 

(44) IU/,||-i/2,c,r = sup —- < sup —- < sup —- = UhW-mx- 

The left inequality is then proved. For the right one, we only have to note that P)j||-i/2,r ^ h^'^(l + 
I log /z|)||/l^Jlo,r' ^rid the inequality follows from the triangle inequality. n 

Lemma 29. It holds for Aj, e S/, that 

(45) ||^||i/2,r < ll^/,lli/2,c,r < (1 + |log/i|)||^||i/2,r 
and 

(46) P/,ll-i/2.r < (1 + \logh\f\\AHU/2,cx- 

Proof For any Ai, e Sh, \\Ah\\i/2,r < llf=i UhXi\\i/2,r = IU/Jli/2,r,r and P/,||i/2,c,r = 2^i P/i;frlli/2,r ^ 
(1 + I log/z|)||/l/,||i/2,r- Thus, (l45l) is proved. 

We now turn to (|46] ). For any ///, e S/,, by Lemma 123) and (|45l) . \\p.h\\i/2,c,r + \\l^l\\i/2,r ^ (1 + 
|log/z|)2||///,||i/2,r- Therefore, 

(Ah,iUh)r < IU/JI-i/2,c,rll)"Alli/2,c,r + ll-^ftlloxIKIIo.r ^ ll^ll-i/2,.-,rll^7lli/2,f,r + /^^^^ll-^ftlloxIKIIiRr 
^ Uh\lu2,cA\\Mi/2,c,r + IKIIi/2,r) ^ (1 + |log/z|)^||/l/,|Ui/2,f,rlNlli/2,r' 

and P/JUi/2,r = ^^P^,i,es„ ^S^ ^ (1 + I log/?|)^PAl|-i/2,c,r- This finishes the proof. n 

5.5. Discrete harmonic extension operator and the discrete inverse generalized Poincare- 
Steklov operator. For A^ e S/„ let Uhj,^ satisfy 

.^j. I (Vmm„,Vv/,) = 0, Vv/, 6y,,o, 

\ Uh,Ai, = A, on5D. 

Define the discrete harmonic operator £■/, : S/, -^ Vhhy EhAh = "Ma- Simultaneously, (fi/j/l/,, A/,jV/,) 
for Vh 6 V/,(). The stability properties of E and E^ can be inherited by E^. The lemma below can 
be referred to [|T3l . 

Lemma 30. ||£'/,/lft||i,n = P/illi/2,r/or^/, 6 S/,. 

Define Si^, := JA G Sh : fi^/, 6 KiA)]. Then, £/,S/, = £,,S/, n ^(A), and by Lemma [H 
codim(!Bh,Sh) = codim(EkBk, EhS/,) = codim{EhSh r\'R(A),Ei,!Bh) < mo. 
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Lemma 31. The following stability results hold: 

(1) it holds for Ah 6 S^ that 

\\EhAii\\[iii > ||/l/,||-i/2,c + h ||/l/,||o,r; 

(2) it holds for Ah 6 S/, that 

\\EHA,Mo,n^Uh\\-i/2,c + h''^\\Ah\\o,r. 
Proof Given Ah e S/„ let v e H^iQ.) D H^q(Q.), then 

dv 

\(,E,Ah - EhAh, Av)| = \(Ah, — )r - (Eh^h, Av)| = HVEhA,,, Vv)| 

= inf KVEhAh, V(v - v,))| < /i||£/,^„||i,n|v|2,n 

< /?'^'||4IMIAv||o,n. 
First, taking v such that Av = E^Ah gives 

II^A^Allcn > CiP,,|Ui/2,, - C2h^'^\\Ah\\o,r, 
and further, H^/i/lAllo.n > c-3/z'^^||/l/,||o,r- It is obtained that 

\\EhAh\\o,n > max(-^i^|U,||_i/2,c.,C3/z^^'lU/,||o,r)- 
C2+C3 

For Ah 6 Sh such that EhAh e "RCA), taking v such that Av = EhAh yields 

l|£A4llo,n< 11^11-1/2,. + /^'^'ll^/Jlo,r- 
The lemma is proved. n 

Define the discrete inverse generalized Poincare-Steklov operator 5/, : S/, -^ S/, by 

(48) Uh,Sh7h)r = (EhAh, Ehjh), V Ah, jh 6 Bh. 

The decomposition of the first biharmonic problem can be inherited in the discrete level. 

Proposition 32. Let Uh 6 V/jo solve (A^^Mft, A/,jV/,) = if,Vh), Vv/, e V/io, ?/zen m/, can be obtained 
by seeking (Uh, ^/,, w/,) e V^o x S/, x y/,o, i'wc/? that 

(49) (AhaUh, \2Ph) = (f, Ph), ypheVho; 

(50) (4,5/,7/,)r = -(A/,,2M/„£ft7/,), 'ijh&Sh, 

(51) (A/,,2M/„ A/,,2Vft) = (A/,,2W/! + £^/i^ft, A/,,2Va), V Vft e V/.o. 
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Define Tu '■ ^h -^ ^h such that for y^,///, 6 Sh, 

(52) (T'h7h,l^h)r ■ = (yEhjh, ^Ei,iuh) + {Ei,y,„ EhUh), 
and Dh : Sh -^ Sh such that for yh,pih 6 Sh 

K 

(53) (2)/,r/„ ///,)r : = J]iVE,0^,r,), VE.O^huh)). 

i=l 

Lemma 33. fh is an isomorphism from (Sh, IMIi/2,r) onto (S/,, IMI-i/2,r), andDh is an isomorphism 
from (Sft, II • ||i/2,c,r) onto (S^, || • ||-i/2,c,r)- 

Some equivalent description of 5/, can be established by means of D/, and Th- 

Theorem 34. For any A^ 6 S/„ it holds that 

(Ah,ShAi,)r > {Dl Ah, /l/,)r + K^l, ^^)r, 

and for any Ah 6 S/,, 

{Ah,ShAh)T < ifDtXr, A„)r + KK, 4)r- 
Proof. The theorem follows from Lemmas |3T] [271 and |33] directly. n 

Theorem 35. For any Ah 6 S/,, it holds that 

(^, 5„^)r > (1 + I log h\r\{Th'Ah, ^/,)r + K^l, ^^r), 

and for any Ah 6 S/,, 

(^,„ S„i„)r < {Th^Ah, A,Or + /?(^^, ^Dr- 
Proof. The theorem follows from Lemmas |29l |271 and |33] directly. n 

Theorem 36. For any Ah 6 S/,, zY /loW* that 

(^,,S,ifi)r > (1 + |log/i|)-4(!r;,-i^,,^,)r, 

and for any Ah 6 S/„ it holds that 

(Ah,S hAh)r ^ ("Fh Ah,Ah)r- 
Proof. The theorem follows from Theorem |35] and Lemma [ 
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5.6. A modified stable decomposition of M/jj. Define V/,o := [ph 6 V^o : (A/,jp/,)|r e S/J and 

Mit,\ := Iwh 6 Mftj : //,w/, e y/,oi. Then, by Lemma[5l codim{Vho, Vho) < codim{!Bh,!Bj,) < mo and 
codim{Mhi, Mhi) = codim{Vho, IhMh^i) < codim(Vho,VhQ) < mo. 

Lemma 37. Le? v^ e y/,o. Let A/,j an J A^ 2 ^^ Laplacian operators defined on Vm- It holds that 

(1) ||A/,,2V/J|o,n < l|A,,,iV/J|o,n; 

(2) \\Ah.iVh\Wii ^ h~^'^\\Ah^2Vh\\o,nforvi, e V^. 

Proof. Let qh 6 Vho, then by the definitions, (A/,jV/,,^/,) = (A/,_2V/„^/0; that is, A/,_2V/j 6 Vho is the 
L^ projection of A/,jV/, into Vho- The first result follows then. 

Now let "H/, be the //'-projection operator from V/, to Vi,o, namely (Vp;„ Vqi,) = (Vpf,, VHhqh), 
for V Pi, 6 Vf,o and 'i qf, £ Vf,. It follows that q,, - 'Hqi, = E,,qi,\Y and (A/,jp/„ ^,,) = (Af,^2Ph, 'J^hQh)- 

Ifqhlr e S/„ then WHhqi, - qhWo.n < ll'?/,lrll-i/2,c,r + /2^''^ll'?/,lrllo,r < ll<?/Jrllo,r < h~^'^-\\qh\\o.n, and thus 
\\H,qh\\o,n < \\qh\\o,n + WHiGh - qhlWii ^ h-^'mqh\\o,a. Therefore for wj, e Vm, (Ah.iWh, Ah^Vh) = 
(A;,,2>v/„'K/,(A/,,iV/,)) < ||'K/,(A/,,iV/0llo,nl|Aft,2W/,||o,n ^ /■?"'^^l|Ai,/,v/,||o,nl|Aft,2W/,||o,n, and ||A/,jW/,||o,n = 

(A/, iW/,, A/, iV/,) _,/j 
sup — -p- -^ <h ' ||A/j2W/J|ofi. The lemma is proved. n 

„„ey„o l|A/,,lV/,||o,n 

Lemma 38. The modified stable decomposition holds: 

(1) For Wh 6 M;, 1, it holds that 

(54) \wh\ljj ^ 2_j ^r^^ll^/' ~ nftj/ftWftllo^r + (Ahahwh, Ah^hwh). 

(2) Forwh 6 Mj,i, it holds that 

(55) \wh\l^h ^ 2-j ^T^^^'^h - T^KihwhWlj + h'\AhaIhWh, Ah^hwh)- 

Ten 

Proof. Combining Theorem [8] and Lemma |37] leads to the lemma directly. n 

6. Optimal solvers for fourth-order finite element problems 

In this section, we will present several eff"ective solvers for both the first and the second bi- 
harmonic problems. We will provide detailed theoretical analysis for these solvers as well as 
numerical experiments that support our theories. 

In the presentation below, Rh : M^' ^ -^ Mh^k represents any symmetric smoother (such as the 
Jacobi and the symmetric Gauss-Seidel smoother) for the discrete biharmonic operators, and Iih,k 
follows from the definition as in Section 1421 
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6.1. Second biharmonic finite element problem A/,2M/, = //, on M/,2. For the second bihar- 
monic finite element problem, the discrete Laplacian operator of second kind provides an optimal 
preconditioned By Theorem [8l Lemma |3] and the FASP theory, we obtain the following theorem. 

Theorem 39. Define 

(56) 5/,2 := Rh + n,,2A;;5(n,,2)*, 

then the m^-th effective condition number /cJ;J|j(5/, 2^/, 2) ^ 1. 

Let Np be the number of interior points of the triangulation, and let Nh2 be the number of the 

degree of freedoms of Ah^2Uh = fh, then A^/,,2 = Np. 

Lemma 40. / [751I47]/ The linear element problem for the Poisson equation can be solved in the 
complexity ofO(Np log Np). 

Theorem 41. When PCG is applied on Bh^^h.i^h = Bh^fh, the total complexity is 0{Nh,2 log A'^;,,2)- 
Proof. The theorem follows from (|9]). Theorem l39l and Lemma l40l n 

6.2. First biharmonic finite element problem A/, jw/, = //, on M/, j. 

6.2.1. Preconditioning effect of the discrete Laplacian operator of the second kind. The discrete 
Laplacian operator of second kind induces a simple preconditioner for first biharmonic finite ele- 
ment problem. 

Theorem 42. Define 

(57) B',y.= R,+n,^,^ll(Ji,^,r. 

then the m^-th effective condition number x^nS^'h i^Ki) ^ h~^. 

Proof The theorem follows from Lemma |3l and Lemma |38l n 

The computational cost is dominantly contained in Poisson solvers. Let A^/, 1 be the number of 
the degree of freedoms of Ah^ui, = fh, then Nh,i = Nh^2 = h~^. 

Theorem 43. The complexity of PCG applied on B'^ jA/j jW/, = 5^ j//, is 0(Nl-^^ \ogNh,i). 

6.2.2. An optimal preconditioner. 
Theorem 44. Define 

(58) Bh,i = Rh + n,,i(A;iA/,i)"'(n/a)*. 
The condition number ofBijj^Ahi is bounded uniformly. 
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Proof. The theorem follows from the FASP theory and the stable decomposition (fT6l) . n 

The main work of the preconditioner 5/,_i is in the inversion of A* ,,Ai./,: namely, given //, e Vm, 
we find Uh e V/,o such that 

(59) (A/,jMft,A/,jv,J = (//„v/,), VvftGyfto- 

Strategy of solving (l59l) . Let w/, be the solution of (|59l) . and vv/, := A/, jM/,, then {uh,Wh) is the 
unique solution of the equation 

(Vm/„ Vp/,) = (wh,Ph) ^ Ph e V/, 



^^^^ ^ (VwA,Vv,) = (/,v,) Vv,ey,o. 

Namely, (|59| ) is equivalent to (l60l) . whereas the latter is the discretization of the mixed formulation 
of the first biharmonic problem given in [jTOl . which finds (m, vv) 6 Hq{Q.) x H^(Q.), such that 



(61) ^ r^- r7„x _./■ „x w„^ r/l 



(Vm, Vv) = (vv, v) Vv6//'(Q), 
(Vw,V;?) = (/,;?) ^pGH'.iQ). 



The coupled system (|6T1) can be decoupled as demonstrated in the following lemma. 

Lemma 45. Let (u, w) be the solution of (|6TI) . then they can be obtained by seeking (u, w. A) G 
//iJ(Q) X H^^{n) X H^'\Y) such that 

(62) (Vw,vp) = if,p), yp€Hl,(ny, 

(63) (EA, Ey) = -(w, Ey), V y 6 //^/'(r); 

(64) iVu,Vv) = {w + EA,v), Vv6//^(Q). 
AnJ, (MjVv) = (m,w + £■/!). 

Proof. Let /I = vv|r, then vv = vv + fi/l, with vv 6 Hq(Q.) uniquely determined by (|62|) . Since 
(vv, £7) = -(Vu,VEy) = for 7 e H^^^iT), (EA,Ey) = -(vv, £7). Thus, A can be solved from 



63I) . Further, w can be solved from (|64l ). The proof is finished. n 



Analogously, the decoupling of (1591) can be carried out by means of E^. Similar to Lemma 
we can prove the lemma below. 

Lemma 46. Let Uh be the solution of (|591) . then it can be obtained by seeking (uh, w/j, A/i) e V/,o x 
V/jo X Sh such that 

(65) {Vwh, Vph) = (/, Ph), V ph e Vi,o; 

(66) (EhAi,, Ehyir) = -(w/,, E,,yi,), V y/, 6 S/,; 

(67) (Vuh, Vv/,) = (w/, + E,,Ah, Vh), V v/, e Vho- 
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According to Lemma |40l (|65] ) and (|67T ) can both be solved optimally. The difficulty lies in 
solving (l68l) in optimal complexity with respect to the size of (|591) . The problem (|66l) is equivalent 
to finding A^ 6 Sh, such that 

(68) iAh,Sh'yk)r = -{wh,Eh-yh), ^yh^^h- 

Optimal and nearly optimal preconditioners for S/,. Based on the equivalent description of 5/i 
in Section 1531 we present three preconditioners for 5/,. They are 

Thy. =IDhI* +h ^JJ*; 

ThX' =iThr+h-'jr; 
Thy- = 'Fh- 

Here, / is the inclusion operator from S/, to S/,, J is the inclusion operator from S^^ to Sh, and /* 
and J* are the adjoint operators of / and J, respectively. 

Theorem 47. It holds for j = 1,2,3 that K^J^^iT^jSh) < 1 + \\og hfj, with pi = 0, jSj = 4 and 
A =4. 

Proof. The theorem follows directly from the FASP theory. Lemma |3] and Theorems l34l |35] and 
[36l respectively. n 

Let Ns be the number of the degree of freedoms of (l66T i. Then A^/, j = A^^, and Nh.i is equivalently 
the number of the degree of freedom of (|59] ). The theorem below follows from Theorem l47l 
Lemma gOl and ^. 



Theorem 48. (1) The equation (|66|) can Z?e solved by PCG with Thi(or T/, 2, T/, 3J a^ f/ie pre- 

conditioner in the complexity ofO{Nh,i \ogNh,\), or 0{Nh,i log^ Nh,\), 0(Nh,i log^ Nh,i), re- 
spectively. 
(2) The problem (l59l) can Z7e solved in the complexity ofO{Nk,i \ogNh,i). 

Evaluating the action of the operator (A^^Ah^y^ requires solving the mixed finite element 
problem (1591 ) once. Therefore the preconditioner 5/, 1 can be carried out in the complexity of 
0{Nh,i log Nh,\). The theorem below is then derived. 

Theorem 49. When PCG is used to solve the problem BhjAh^iUh = Bh^ifh, the total complexity for 
convergence is 0{Nh,i \ogNh,\). 

7. Numerical examples 

In this section, we present several numerical examples to illustrate the eff'ects of the precondi- 
tioners given in last section. We will compute and record the extremal eigenvalues of the precon- 
ditioned operators. We will also test the performance of PCG method with given preconditioners 
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on some modal problems. We run the various PCG computations with the starting guess 0, and 
with a stop criteria whereby the relative residual (||residual||/2/||rhs||/2) is smaller than 10"^. 

We test the preconditioner on both convex and nonconvex domains, as shown in Figure [U We 
divide each computational domain by successively refined quasi-uniform meshes, and we carry 
out numerical experiments on the multiple meshes to test the performance of each preconditioner. 
In the tables below, we use A for an eigenvalue, k for a condition number, and DOF for the number 
of the degree of freedom. 
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Figure 1 . The computational domains with initial mesh. 

7.1. Preconditioning effect of 5/, i for first biharmonic finite element problem. We take the 
frequently used Morley element ( [|24]| ) and the modified Zienkiewicz element ( ll36l ) as examples 
to demonstrate the preconditioning effect of the preconditioner 5/, j with respect to A/, j. We use 
the symmetric Gauss-Seidel method with three iterations as a smoother. 

As shown in Tables [H [21 [H and|4l we observe that for both the Morley element and the modified 
Zienkiewicz element and on both convex and nonconvex computational domains, the condition 
number of 5/, jA/, j is uniformly bounded as the triangulation is refined. 



nodes 


Morley Element 


Modified Zienkiewicz Element 


DOF 


^1 


^N,a 


K 


DOF 


^1 


^W/,.i 


K 


961 


3969 


0.70 


1.61 


1.1% 


2883 


0.18 


1.82 


10.08 


3969 


16129 


0.70 


1.61 


2.30 


11907 


0.18 


1.82 


10.06 


16129 


65025 


0.70 


1.61 


2.29 


48387 


0.18 


1.82 


10.12 



Table 1. Eigenvalues of 5/j_i A/, i: unit square. 



nodes 


Morley Element 


Modified Zienkiewicz Element 


DOF 


^1 


M.i 


K 


DOF 


^1 


^N,,,x 


K 


2089 


8529 


0.70 


1.61 


2.30 


6267 


0.16 


1.74 


10.68 


8529 


34465 


0.70 


1.60 


2.29 


25587 


0.16 


1.74 


10.83 


34465 


138561 


0.70 


1.61 


2.29 


103395 


0.16 


1.74 


10.90 



Table 2. Eigenvalues of 5^1 A;, i : hexagon. 
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nodes 


Moiiey Element 


Modified Zienkiewicz Element 


DOF 


Ai 


Ma 


K 


DOF 


Ai 


^Nia 


K 


3201 


12033 


0.57 


1.61 


2.83 


8835 


0.16 


1.68 


10.55 


12545 


48641 


0.56 


1.61 


2.86 


36099 


0.16 


1.68 


10.53 


49665 


195585 


0.56 


1.61 


2.88 


145923 


0.16 


1.69 


10.55 



Table 3. Eigenvalues of 5/jiA/ji : "L-" shape domain. 



nodes 


Morley Element 


Modified Zienkiewicz Element 


DOF 


Ai 


Anu 


K 


DOF 


^1 


An,,,, 


K 


2673 


10017 


0.70 


1.61 


2.29 


7347 


0.16 


1.65 


10.31 


10465 


40513 


0.70 


1.61 


2.30 


30051 


0.16 


1.65 


10.34 


41409 


162945 


0.70 


1.61 


2.30 


1539 


0.16 


1.66 


10.35 



Table 4. Eigenvalues of 5/, i A/, i : trident domain. 

7.2. Preconditioning effect of 3^2 for the second biharmonic finite element problem. We take 
the Moriey element ( (231, Il35l ) to demonstrate the preconditioning effect of the preconditioner 
5/, 2 with respect to A/, 2- Numerical verifications ( [|35l ) showed that the Morley element fits the 
second biharmonic problem. 

As shown in Tables [51 [6l and|71 we observe that for both convex and nonconvex cases the num- 
ber of PCG iterations remains nearly constant as the number of DOF grows. Moreover, when 
the computational domain is convex, the eigenvalues of the preconditioned operator are bounded 
uniformly. However, when the computational domain is not convex, the eigenvalues of the precon- 
ditioned operator are bounded uniformly from below. In addition, when the domain is not convex, 
the eigenvalues are bounded uniformly from above with no more than mq exceptions. 



Refine 
Times 


Four-Sc 


uare 


Hexagon 


"L-"shape domain 


Trident 


DOF 


step 


DOF 


step 


DOF 


step 


DOF 


step 


5 


16385 


14 


34817 


11 


49153 


17 


40961 


12 


6 


65537 


14 


139265 


11 


196609 


17 


163841 


12 


7 


262145 


15 


557057 


10 


786433 


19 


655361 


12 



Table 5. PCG steps needed for solving 5/,,2A/,,2- 



Refine 
Times 


Four-Sc 


uare 


Hexagon 


DOF 


^1 


Ani,^ 


K 


DOF 


^1 


Ani,i 


K 


5 


16385 


0.71 


1.60 


2.25 


34817 


0.75 


1.57 


2.09 


6 


65537 


0.71 


1.60 


2.25 


139265 


0.78 


1.53 


1.96 


7 


262145 


0.71 


1.60 


2.25 


557057 


0.81 


1.51 


1.86 



Table 6. Eigenvalues of 5;,2A/,2: convex domain. 
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"L-"shape Domain(mo = 1) 


Trident Domain(mo 


= 2) 


DOF 


^1 


^Ns-l 


An^ 




DOF 


^inin 


Ans-2 


A-Ns-l 


Ans 




193 


0.73 


1.61 


1.87 


2.22 


161 


0.69 


1.61 


1.64 


1.68 


2.34 


769 


0.69 


1.61 


2.33 


2.33 


641 


0.70 


1.61 


1.67 


1.72 


2.29 


3073 


0.71 


1.62 


3.17 


2.28 


2561 


0.71 


1.62 


1.80 


1.88 


2.28 



Table 7. Eigenvalues of Bh^^Ki- nonconvex domain. 

7.3. Preconditioning effects of the preconditioners for 5/,. 

The optimal preconditioner T/, j. We record the performance of T/, j for S/, in Tableland Table 
[9l We observe that, the eigenvalues of Th,\Sh are uniformly bounded from below. In addition, they 
are bounded from above with no more than mq exceptions. 



Unit Squai-e 


Convex Hexagon 


DOF 


^1 


"^Ns 


K 


DOF 


^1 


'^Ns 


K 


256 


0.20 


4.48 


21.96 


352 


0.21 


4.87 


23.40 


512 


0.20 


4.49 


22.52 


704 


0.20 


4.89 


24.13 


1024 


0.20 


4.49 


22.86 


1408 


0.20 


4.90 


24.62 



Table 8. Eigenvalues of T/,! 5/,: convex domain. 



"L-"shape Domain(mo - 1) 


Trident Domain(m() - 2) 


DOF 


^1 


Ans-I 


An^ 


^ett 

""JO 


DOF 


"-min 


Ans-2 


Ans~\ 


M- 


ett 


512 


0.21 


4.50 


35.37 


21.89 


448 


0.19 


4.50 


14.70 


17.08 


24.07 


1024 


0.20 


4.50 


56.70 


22.39 


896 


0.18 


4.50 


20.24 


23.39 


24.63 


2048 


0.20 


4.50 


90.58 


22.76 


1792 


0.18 


4.50 


27.59 


31.75 


25.01 



Table 9. Eigenvalues of r/,iS/,: nonconvex domain. 



The nearly optimal preconditioners Th2 and T/, 3. We also carry out the same numerical experi- 
ments for the performance of T/, 2 and T/j 3 as preconditioners for Sh, recording the results in Tables 
llOll2llll and[T3l We observe that the eigenvalues of Th,2Sh and T/i^S/, are bounded from above, 
with no more than mo exceptions. In addition, they are bounded from below, with slight decrease. 
To illustrate the difference and make a comparison, we compute the eigenvalues of T^jS /, and 
Th,T,S k on all the computational domains and plot the distribution of them in Figure |2] We observe 
in Figure |2] that, T/, 2 performs better than T/, 3 at capturing the extremely low-frequency part, and 
is as good as 7/^3 at capturing the high-frequency part. 

8. Concluding remarks 

This paper is aimed at developing efficient iterative methods for solving the algebraic systems 
arising from direct finite element discretization of the boundary value problems of fourth-order 
equations on an unstructured grid. The following objectives have been accomplished: 
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DOF 


Unit Square 


Hexagon 


M 


Ans 


K 


Ai 


Ans 


K 


512 


1.22 


9.34 


IM 


0.29 


4.88 


16.63 


1024 


1.07 


9.34 


8.72 


0.25 


4.90 


19.73 


2048 


0.95 


9.34 


9.87 


0.21 


4.91 


23.38 



Table 10. Extremal eigenvalues of Tu^Sh- convex domains. 



DOF 


"L-"shape Domain(mo - 1) 


Trident(mo - 2) 


^1 


Ans-\ 


Af^^ 


„ett 


^1 


^^5-2 


Ans-\ 


Aj^g 


^ett 


512 


0.33 


3.25 


65.57 


9.75 


0.19 


3.94 


22.07 


24.66 


20.33 


1024 


0.31 


3.37 


102.92 


10.98 


0.16 


4.01 


29.01 


32.52 


25.24 


2048 


0.28 


3.49 


162.00 


12.35 


0.13 


4.05 


38.07 


42.79 


30.62 



Table 11. Extremal eigenvalues of T/jiS/,: nonconvex domains. 



DOF 


Unit Square 


Hexagon 


^1 


Ans 


K 


^1 


Ans 


K 


352 


0.23 


9.33 


40.64 


0.027 


4.88 


180 


704 


0.21 


9.33 


45.46 


0.023 


4.90 


216 


1408 


0.18 


9.33 


50.48 


0.020 


4.91 


252 



Table 12. Extremal eigenvalues of Th,-iSh- convex domains. 



DOF 


"L-"shape Domain(mo = 1) 


Trident(mo = 2) 


^1 


Ans-\ 


Ans 


ett 


^1 


Ans-2 


Ans-\ 


An^ 


ett 
"mo 


448 


0.067 


3.25 


35.19 


48.53 


0.035 


3.94 


10.20 


11.00 


114 


896 


0.059 


3.37 


54.08 


57.04 


0.030 


4.01 


12.71 


13.83 


136 


1792 


0.053 


3.49 


83.84 


66.14 


0.025 


4.04 


15.93 


17.45 


159 



Table 13. Extremal eigenvalues of and T^^jS i,: nonconvex domains. 



- A class of nearly optimal iterative methods are developed for discrete fourth-order problem 
with both the first and second kinds boundary value conditions in two dimensions. 

- A complete and rigorous analysis is provided for all the algorithms proposed in the paper. 

- Numerical experiments are carried out to confirm all the theoretical results in the paper. 

The algorithms and theories in this paper are valid for general unstructured grids in general polyg- 
onal domains which can be both convex and nonconvex. The iterative algorithms developed are 
the first and the only known methods in the literature for fourth-order problems that are provably 
(nearly) optimal. 

To accomplish the objectives, the Fast Auxiliary Space Preconditioning (FASP) method ( BOl 
|4T]| ) is used as the technical framework for designing the pre conditioners, and the solution of 
a fourth-order problem is reduced to several second-order problems on the discrete level together 
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Figure 2. The distributions of eigenvalues of Tu^Sh and Th^^Su- In each subfigure, 
we use ^ for the eigenvalue of T^^^^, and o for that of Th^-^Sh- 



with local relaxation methods. A number of intricate Sobolev spaces (such as the normal derivative 
trace space of H^(Q.) n H^(Q.)) defined on the boundary of a polygonal domains are carefully 
studied and thereafter used in the analysis of preconditioners. A special observation can be made 
that a straightforward mixed finite element discretization is used as the major component in the 
proposed preconditioners. Indeed, the aforementioned mixed method is either non-optimal or 
non-convergent method as a discretization method ( [[T6l 12111251 [321 I^TlD for the original fourth- 
order problem, and it is interesting to notice that the mixed method provides a nearly optimal 
preconditioner when it is used in conjuncture with additional local smoothers and preconditioned 
conjugate gradient methods. 

Finally, the algorithms and theories in the paper need to be extended to the following cases in 
the future works: 

- General shape regular unstructured grids that are not assumed to be quasi-uniform. 

- Three dimensional case. 

- More complicated fourth-order equations such as Cahn-Hilliard equations. 



Bibliographic comments. The analysis in this paper consists of numerous technical results. Here 
we give a brief description of how some of these results are related to existing results in the 
literature. 

The right inequality of (fT3l ) for first biharmonic problem was first studied in Babuska, Osbom, 
and Pitkaranta [|2l for convex domain in their analysis for some mixed methods, and Hanisch [fT9ll 
for nonconvex domain also for analysis of mixed methods. In the present paper, we establish (fUl l 
for both first and second biharmonic problems on both convex and nonconvex domains. 

Theorem [T5l can also be found in different form in Peisker [|27l . who gave a proof assuming 
the domain is convex. In the present paper, we establish Theorem [15] and prove it by means 
of an auxiliary Stokes problem on both convex and nonconvex domains. Lemma [20] presents 
the isomorphisms between the normal derivative trace space and the Laplacian trace space of 
biharmonic functions on polygonal domains. When Q. is smooth, similar result was given by 
Glowinski and Pironneau [[T3l. 
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Peisker [|27l studied Lemma [311 and Braess and Peisker [|3l proved Lemma [37] both in a special 
case that Q is convex and by analyzing the property of the discrete extension operator E^ associ- 
ated with the harmonic extension operator E. In the present paper, we study the property of Eh 
associated with the generalized harmonic extension operator E^, which does not coincide with E 

in H'^^^^iY) \ H^l^iX) or in nonconvex domains, and establish Lemmas [3T[ and [TT] for both convex 
and nonconvex domains. A disguised form of Lemma [46] can also be found in Glowinski and 
Pironneau [[T3]| . but our formulation and proof are quite different. 

Our preconditioner ([58] ) for the first biharmonic problem was motivated by an algorithm pro- 
posed in Peisker and Braess [[281 . where an algebraic preconditioner for Morley element problem 
was presented for the special case that the domain is convex and the triangulation 7"/, consists of 
triangles that are similar to each other. Peisker [[271 also noticed the role of 5^ and presented a 
preconditioner for S u in matrix form when Q. is convex. Her preconditioner can be realized by 
fast Fourier transform (FFT) on a graded bisection mesh of dQ.. The technique of analyzing an 
interface operator by Fourier analysis was also used in, e.g., [[9] [HI. In the present paper, we 
design preconditioners in a unified framework for both first and second biharmonic problems dis- 
cretized by various finite element methods on general triangulation for both convex and nonconvex 
domains. And, we establish preconditioners for 5/, in a more analytic approach, for general polyg- 
onal domain Q that are triangulated by general meshes. 
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